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^ , Abstract 

In this paper we derive expressions for matrix elements (4n,H<j)j) for the Hamiltonian H = 
— A + ^2 q a(q)r q in d > 2 dimensions. The basis functions in each angular momentum subspace 
are of the form <fii{r) = r l + 1 +( t - d )/ 2 e -" rP / 2 j i > 0, p > 0,i > 0. The matrix elements are given in 
terms of the Gamma function for all d . The significance of the parameters t and p and scale s are 
discussed. Applications to a variety of potentials are presented, including potentials with singular 
repulsive terms of the form f3/r a , a,f3 > 0, perturbed Coulomb potentials —D/r + Br + Ar 2 , 
and potentials with weak repulsive terms, such as — ^r 2 + r 4 , 7 > 0. 

PACS 03.65. Ge, 31.15.Bs, 02.30.Mv. 
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1. Introduction 

We study quantum mechanical Hamiltonians H = —A + V(r) in d > 2 dimensions, where 
V is a spherically-symmetric potential that supports discrete eigenvalues, r = |r|, re 3? d . We 
estimate the spectrum of H in an n -dimensional trial space lying inside an angular-momentum 
subspace labelled by I and spanned by radial functions with the form 

n-1 

4>(r) = Y,<V- i+1+(t - d),2 e- hr \ *>0. (1.1) 

i=0 

If the potential is chosen to be a linear combination of powers 

V(r) = 5>(«)r«, (1.2) 

then all the matrix elements of H may be expressed explicitly in terms of the Gamma function. 
The expressions obtained will be functions of the parameters t and p , and also of a scale parameter 
s to be introduced later. If the potential V is highly singular, the parameter t must be chosen 
sufficiently large so that (V) exists. The advantages of the particular form chosen for the radial 
functions will become clear in the development. Thus we have n + 3 variational parameters with 
which to optimize upper estimates to the spectrum of H , with one degree of freedom being employed 
for normalization. 

Systems with Hamiltonians of this type have enjoyed wide attention in the literature of quantum 
mechanics [1-62]. This interest arises particularly from the usefullness of these problems as models 
in atomic and molecular physics. Many numerical and analytical techniques have been used to 
tackle Hamiltonians of this form. In Section 2 we derive general matrix elements and show how the 
minimization with respect to scale s can be easily included. In section 3 we discuss some numerical 
issues not the least of which is the usefulness of the reduction of the matrix cigen equations to 
symmetric form by first diagonalizing the 'normalization' matrix N — \{4>u4>j)]- The dependence 
of the eigenvalues on the parameters {p, t, s} may be rather complicated. Since changes to scale s 
do not involve the recomputation of the basic matrix elements, a policy which emerges is to fix n, 
always optimize fully with respect to scale s, and, if necessary, optimize approximately with respect 
to t and p by exploring a few values; if higher accuracy is required, a full optimization is undertaken, 
or n is increased. In Section 4 the matrix elements are applied to a variety of problems and the 
results are compared with those found in earlier work. We suppose that the Hamiltonian operators 
in this paper have domains T>(H) C L 2 (3? d ), they are bounded below, essentially self adjoint, and 
have at least one discrete eigenvalue at the bottom of the spectrum. This, of course, implies that the 
potential cannot be dominated by repulsive terms. Because the potentials are spherically symmetric, 
the discrete eigenvalues can be labelled by two quantum numbers, the total angular momentum 



A basis for variational calculations in d dimensions 



page 3 



£ = 0, 1, 2, . . . , and a 'radial' quantum number, n = 0, 1, 2, . . . , which counts the eigenvalues in each 
angular-momentum subspace. These eigenvalues satisfy the relation E^ e < E^, n < m. With our 
labelling convention, the eigenvalue E d ni (q) in d > 2 spatial dimensions has degeneracy 1 for £ = 
and, for I > 0, the degeneracy is given [63] by the function A(d,£) , where 

A(d,£) = (2£ + d-2)(£ + d-3)\/{£l(d-2)l}, d>2,£>0. (1.3) 

Many techniques have been applied to approximate the spectrum of singular potentials of the form 
(1.2) using perturbation, variational, and geometrical approximation techniques [1-62]. Exact solu- 
tions for the energy may be obtained in some special cases by first choosing a wave function with 
parameters, and then finding a potential of the form (1.2) for which this wave function is an eigen- 
function; this is possible only when certain constraints are satisfied between the parameters {a(q)} , 
as we shall discuss later. 

2. Matrix elements 

We consider first the action of the Laplacian in d dimensions on a wave function *(r) = 
tp(r)Ye(8o, #i, • • • , with a spherically-symmetric factor ip(r) and a generalized spherical har- 

monic factor Ye,. If we remove the spherical harmonic factor after the action of the Laplacian on $ 
we obtain [64] 

A * ,/// n d-1 ,,. . £(£ + d-2) ,. . . ,. 

Ye r r z 

The radial Schrodinger equation for a spherically symmetric potential V(r) in d -dimensional space 

is therefore given by 

~W ^ t + ~ ^ ^ + V ^ - E ^ ^(r)eL a ([0 l oo) > r-- 1 dr) (2.2) 

A correspondence to a problem on the half line in one dimension with a Dirichlet boundary condition 
at r — is obtained with the aid of a radial wave-function R(r) defined by 

R(r) = r^-^VM. d > 2 > R(0) = 0. (2.3) 

If we now re-write (2.2) in terms of this new radial function, we obtain the following Schrodinger 
equation for a problem on the half line 

]2r> 

HR= -— + UR = ER, R e L 2 ([0,oo),dr), (2.4) 
where the effective potential U(r) is given by 

UM-VM+ y + '-W + '-O) , ,2,5) 
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and H = —-§pz + U is the effective Hamiltonian. We note that in (2.5), d and I enter into the 
equation only in the combination 21 + d in U (r) : consequently, the solutions for a given central 
potential V(r) are the same provided d + 21 remains unaltered. In this setting, our trial wave 
functions now have the explicit form 

R,(r) = r {t+1)/2+l exp(-r p /2) G L 2 ([0, oo), dr). (2.6) 

Thus we have for the general radial function in our trial space 

d-l 
i=0 

The matrix elements we seek (in a given angular momentum subspace) are given by 

Ha = (Ri, -BfJ) + J2 a (l) ( R i, r9R i) ■ ( 2 - 8 ) 
<? 

For each potential term r 9 , if we everywhere omit the constant angular factor (equal to At: in the 
case d — 3), we find the following fomulae, expressed now in terms of the £ 2 ([0, oo), dr) inner 
product: 

/■OO 

P ij (q,p,t) = (R i ,r^R j )= r'+W+^dr, i,j = 0, 1, 2, . . . , 

(2 9) 



P V P 

This type of integral is found by setting x — r p , and using the differential relation r k dr = 
(l/p)x( k+1 ~ p ^ p dx and the definition of the Gamma function. The normalization integrals are special 
cases of (2.9), namely 

N l3 (p,t) = (Ri,Rj) = P«(0,p,t) = *r ( l + J+ J + 2 ) ■ (2.10) 



P \ P 

After some algebraic simplifications we find that the corresponding kinetic energy matrix elements 
Kij (p, t) = - (Ri,R") are given by 

K aM = ^ T { l + j p + t ) [( 2 £ + d - 1 )( 2i + d -3) + i-(i-3) 2 +p(i + j + t)] , t>0. 

(2.11) 

We note that these terms of the Hamiltonian matrix elements Hij are all symmetric under the 
permutation (ij) (because of Hermiticity), and invariant with respect to changes in d and t that 
leave the form 2£ + d invariant. These formulae may be used as they stand for all dimensions d > 2 
provided that t > is chosen sufficiently large t > — (2 + q) to control the most singular potential 
term . We note, in addition, that the choice {d = 3, I = 0} also provides the odd-parity solutions 
in one dimension. 



A basis for variational calculations in d dimensions 



page 5 



We now consider the problem of minimizing (R, HR) with respect to the vector v of coefficient 
{ci}™~Q subject to the constraint that (R, R) — 1 . We immediately obtain the necessary condition: 

Hv = £Nv (2.12) 

By the min-max characterization of the spectrum [67], the eigenvalues of this matrix equation are 
upper bounds to the unknown exact eigenvalues Eu,i = 0, 1,2, .. . n — 1. We assume that these 
discrete eigenvalues of the underlying operator H are either known to exist, or indeed are demon- 
strated to exist by the results of this variational estimate. By considering scaled radial wave functions 
of the form 

R s (r) = R(r/s), (2.13) 

we find that factors of s remain only according to the dimensions of the terms. In effect, when using 
the scaled wave functions (2.13), we can leave the matrix N unchanged and replace the matrix for 
H by 

Hij(s) = ^Kij( P ,t) + J2a(q)s q P ll (q,p,t). (2.14) 
Thus the upper bounds we seek are provided by the eigenvalues of the matrix equation 

H(s)v = £Nv, (2.15) 
which now depend, for a given n and £, on s,p and t and we write 

E u <£ u = £ u (p,t,s), i = 0,l,2...n-l. (2.17) 

The problem now is to find these upper estimates and minimize them with respect to the three 
parameters {p, t, s} . 

3. Some numerical considerations 

Rather than solving the general matrix eigenequation (2.7) directly, it is often desirable to 
use the fact that N is positive definite to transform the problem to symmetric form. In physics 
literature this is sometimes called a Lowdin transformation [65] and is equivalent analytically to 
converting the basis functions to an orthonormal set by applying the Gram-Schmidt procedure. We 
first diagonalize TV with the aid of an orthogonal matrix, say S. We then get S T NS = M~ 2 : the 
square root M exists because N positive definite, which implies that the diagonal matrix has only 
positive eigenvalues. The original problem (2.12) (or the scaled version (2.15)) may now be written 
as 

Hv = \Nv -» S T HSS T v = \S T NSS T v = XM- 2 S T v (3.1) 
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If we multiply on the left by symmetric diagonal matrix M we obtain 



MS T HSMM~ 1 S T v = XM~ 1 S T v 



(3.2) 



If we now write H — MS T HSM 1 and u — M 1 S T v, we obtain the reduction 



Jiu = Am, 



(3.3) 



where H T — H. This is the symmetric alternative to our original eigenvalue problem. It has also 
been shown that the Cholesky decomposition [66] in which the matrix N is written N = L T L, 
where L is upper triangular, is often numerically faster and more stable than finding the square 
root M . Computer algebra systems often allow one to solve these problems directly without knowing 
which method is in fact implemented; the main purpose of our remarks is to show constructively 
that solutions are always possible. 

Another issue is to do with the Gamma function generating large numbers before (or without) 
the symmctrization of H . To deal with this problem we have found it useful at an early stage to 
divide all the matrix elements by (NuNjj)i . 

Ideally the matrix eigenvalues should simply be optimized with respect to the parameters 
{p, t, s}. In practice this is not always a trivially easy task. Typically one chooses the basis di- 
mension n and the angular momentum £, and then finds the n eigenvalues. These numbers must 
be sorted to find, say, the fcth eigenvalue £ki(p,t,s), and finally this function must be optimized 
with respect to the three parameters. This appears to be straightforward until one realizes that the 
matrix eigenvalue problem must be re-solved for each choice of the parameters and, of course, the 
original ordering can be upset. Logically the fcth always has the same numerical meaning but the 
effect is to make the function £ke(p 7 t, s) complicated. It is helpful to note that the basic matrices 
N(p,t), P(q,p,t) and K(p,t) do not depend on s : the Hamiltonian matrix H depends on s by 
the scaling equation (2.10). In order to reduce the difficulty of the search for a minimum we have 
sometimes found it useful to fix p and t and to minimize at first only with respect to s; if neces- 
sary a graph can be plotted of the dish-shaped function £u{s) to give a picture of the minimum. 
This task may then be repeated for some other choices of p and t . In many cases an algorithm 
such as Nelder-Mead tackles the full minimization problem very effectively and there is no more ado 
concerning it. We shall make some comments concerning these matters along with the applications 
described in section 4 below. 

4. Applications 
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We may immediately employ the matrix elements found to solve the eigenvalue problems for 
the general family of Hamiltonians given by 

H = -A + J2a(q)r", (4.1) 
q 

One family we shall study in particular is the class of anharmonic singular Hamiltonians 

^=-^ + - 2 + E^ 'E[0,oo) ( 4 - 2 ) 

q=0 

where a q and X q are positive real numbers, and we assume that the exact wave function tp 
of H satisfies a Dirichlet boundary condition, namely tp(0) = . Inverse power-law potentials 

N 

V(r) = "^2 X q /r aq appear in many areas of physics and for this reason have been widely investigated. 
The spiked harmonic oscillator Hamiltonian, for example, 

H ( a > X ) = -^+r 2 + ^> «>0,A>0 (4.3) 

has been the subject of many mathematical studies which have greatly improved the understanding 
of singular perturbation theory [1,47]. Many different methods [1-62] have been used to study the 
anharmonic singular Hamiltonians (4.2), such as numerical integration of the differential equation, 
perturbative schemes specifically developed for this class of Hamiltonian, and variational methods. 
Among the various methods, the variational method is widely used for calculating energies and 
wave functions since it has the advantage that the eigenvalue approximations are upper bounds [67] . 
Many variational techniques used in the literature were design to solve specific classes of Hamiltonian 
such as (4.3). Aguilera-Navarro et al [31], for example, reported a variational study for the ground- 
state energy of the spiked harmonic oscillator (4.3) valid only for a < 3 . Their study makes use 
of the function space spanned by the exact solutions of the Schrodinger equation for the linear 
harmonic oscillator Hamiltonian, supplemented by a Dirichlet boundary condition ^(0) = , namely, 
ip n (r) = A n e~ r2 / 2 H 2n+ i{r), A~ 2 = 4™(2n+ 1)!^, n = 0,1,2, .. . where H 2n +i(r) arc the Hcrmitc 
polynomials of odd degree. The matrix elements of the operator r~ a ,a < 3 , in this orthonormal 
basis were given as 

-a = , +n V(2m+l)!(2n+l)! r(f) " fc ( m \ T(k + ^)I> + f - k) 
mn 1 ' 2™+«m!n! r(n+f)^ ' \k ) T(k + §)r(§ - k) ' 

A variational analysis was carried out and the ground-state upper bounds were reported for the 
case of a = 5/2. Fernandez [53], soon afterwards, design a particular trial function ip{r) = 
r fc+i e -isr 2 -| tI,2 jS > o, t > 0, to study the ground-state energy of (4.3) for a even integer and 
for arbitrary value of A > . An upper bound to the ground-state of (4.2) was found by a mini- 
mization with respect to {s, t} of E (s, i) = ((1 — t 2 )!^ + 3tl 2 + 2stI- 2 — s 2 I-a + XI 2 - a )/I 2 where 
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I n (s, t) — / °° r" exp(— s/r 2 — tr 2 )dt, — oo < n < oo . These variational results however were not very 
accurate, even for arbitrary large value of A owing to the accumulated error in the computation of 
I n (s,t) . An interesting consequence of Fernandez's work was, however, the exact solution of very 
particular class of (4.3), namely H = —d 2 /dr 2 + r 2 + 9/64r~ 6 , where the exact wavefunction in this 
case reads ip(r) — r%e~ ts t 2 ~i r2 and the exact ground-state energy is Eq = 4 . Aguilera-Navarro 
et al [34] afterwards designed another trial function particularly devoted to analyze the ground-state 
energy of the Hamiltonian H(A, A) . Non-orthogonal basis set of trial wave functions were introduced 
by means of tpn(r) = A n i-Fi(— n; §; r 2 ) exp(— ar 2 — £), n = 0,1,2, ... where A n is the normal- 
ization constant and \Fi (— n; 3/2; r 2 ) is the confluent hypergeometric function. The expressions for 
the matrix elements H mn {&, A) were given by 

H mn = ]T E ^"w^ t( 4m + 3 ^ 2 P + 2 « + 4) + 2 V\I(2p + 2q+ 1) - 4 9 VA(2p + 2« - 1)] 

where the definite integrals I(u) — J °° r u exp(— r 2 — {2^/~\/r))dr were computed by means of the 
recursive relations (u + l)I(u) = (u — l)I(u — 2) + 2\/A7(u — 3) . The shifted factorial (a)„ is defined 

by 

(a) = l, (a) n = a(a+l)(a + 2)...(a + n-l), for n = 1, 2, 3, . . . , (4.4) 

which may be expressed in terms of the Gamma function by (a)„ = T(a + n)/r(a), when a is 
not a negative integer — m , and, in these exceptional cases, (— m)„ = if n > m and otherwise 
(— m) n = (— l) n m!/(m — n)!. The ground-state of (4.3) with a = 4 then follows by diagonalization 
of H in the nonorthogonal basis. This particular study was then extended [36] to provide a global 
analysis of the ground and excited states for the successive values of the orbital angular momentum 
of the super-singular plus quadratic potential r 2 + A/r 4 . Another variational study of the ground 
state of (4.2) was introduced by Hall et al [38] where three parameters trial functions ip(r) = 
r p+e exp(— (3r q ), p = (a — l)/2 were used to approximate upper bounds of the ground-state of (4.3) 
for arbitrary a and A through the minimization of the right-side of the inequality Eq < Eq , where 



EX = min 

e,/3,q>0 



(2p+q+2e-l) gi --(p+e)(p+e-l)g 2 -^g 3 
q Z 



+ 



1 \ 2,q 

-j g 4 +\(2p) a/q g 5 



and 



r (2p+2e+3\ T ( 2 P+ 2e- a+ l \ p ( 2p+2e+l 



U>= y \q y ^ = ^\^— q y 96 -y g 

In attempt to provide a comprehensive variational treatment of the spiked harmonic oscillator Hamil- 
tonian (4.3), for ground-state energy as well for excited states, independent of particular choices of 
the parameters a and A , Hall et al [40-48] based their variational analysis of the singular Hamilto- 
nian (4.1) on an exact soluble model which itself has a singular potential term. They have suggested 
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and used trial wave functions constructed by means of the superposition of the orthonormal functions 
of the exact solutions of the Gol'dman and Krivchenkov Hamiltonian 

d 2 2 A 



H o = -^+r Z + ~- (4-5) 



r 



The Hamiltonian is the generalization of the familiar harmonic oscillator in 3-dimcnsion —d 2 /dr 2 + 
r 2 + 1(1 + l)/r 2 where the generalization lies in the parameter A ranging over [0, oo) instead of 
only values determined by the angular momentum quantum numbers I = 0, 1, 2, ... . The energy 
spectrum of the Schrodinger Hamiltonian H is given, in terms of parameter A as 

= 2(2n + 7 ), n = 0,l,2,..., (4.6) 

in which 7 = 1 + A + j and the normalized wavefunctions are 



Mr) = (-lrW^^rT-ie-i^iFiC-njrr 2 ). (4.7) 



Here i-Fi is the confluent hypergeometric function 

—n)k' 
(b) k k 



iFi(— n; b; z) = , ( n -degree polynomial in z ). (4.8) 



fe=0 

Explicit matrix elements of the Hamiltonian (4.3) can often be found in this orthonormal basis. For 
instance, the matrix elements of the singular operator Ar~ a assume the form 



^ ' h)n Th) V n\m\ 3 2 



1 , (4.9) 



(7)n r( 7 ) V n\m\ '\ 1 ,l-%-n 
where the hypergeometric function 3 F 2 is defined by 

„ / -m,a,b A \^ {-m) k (a) k {b)k , , , . n 

3 ^ 2 ~a 1 = V i / ,V , m - degree polynomial . 

Upper bounds to the energy levels of the Hamiltonian (4.3) then follow by diagonalization of H in 
the orthonormal basis (4.5). In the case where a is a non-negative even number a = 2, 4, 6, ... , the 
hypergeometric function 3 F 2 in (4.9) can be regarded as a polynomial of degree f — 1 instead of an 
m -degree polynomial. Consequently the matrix elements assumes much simpler expressions which 
are useful in numerical computational. For a ^ 2, 4, 6, ... , the variational computational were then 
based on direct use of the matrix elements in terms of the hypergeometric function 3 F 2 . According 
to our discussion up to this point, it is clear that most of the variational methods developed in the 
literature were specifically design to solve the eigenvalue problem of different classes of the singular 
Hamiltonian (4.2). No basis set or trial wave function were design to treat a problem such as the 
singular potentials which at the same time can be used, say, for Hamiltonians with polynomial type 
potentials. The purpose of our basis introduced in section (2) and (3) is to have avaliable at our 
disposal a working variational approach that can be used without a particular references to specific 
potentials or special values for the parameters involve. In the next we apply the matrix elements 
discussed in section (2) and (3) to solve a number of different eigenvalue problems. 
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4.1 Spiked Harmonic Oscillators 

We start our applications by investigating the energy levels of the spiked Harmonic Oscillator 
Hamiltonian (4.3). As we mentioned in section 3, the problem of finding the eigenvalues reduces to 
diagnalizing the real symmetric matrix H = MS T HSM . For a = 2 , the Hamiltonian (4.3) admits 
an exact solutions (4.6). Thus it serves as a benchmark for our variational approach. In Table (1), 
we report our upper bounds for the ground-state of the spiked harmonic oscillator Hamiltonian (4.3) 
for several values of the parameters A and a = |, 1, |, 2, | along with some results obtained in the 
literature. For a = 2 , with m = n = , Table 1 shows that minimization over the three variables 
{p, t, s} yields excellent agreement with the exact solutions (4.6). Such results can be explained by 
observing that direct substitution of the trial wave function ipo(r) into the eigenvalue problem 

H ^ = "^r + ^ + = E o^Po (4.10) 



yields for r — > that 



l — l + ^_(i±^_A = (4ii) 

4 r 2 T 4 r 2-2p 2r 2 ~P r a v^-^-v 



Consequently, for a = 2 and p > , the value t = 1 + 1 1 — \f\ + 4A| yields the best possible 
value of t . As for a < 2 , similar reasoning yields for r — > that t = 1 is an excellent initial 
approximation for t, that is to say, suitable for starting the minimization process. In Table 2, we 
present a comparison between different variational approaches for computing upper bounds to the 
ground-state of the Hamiltonian H = —d 2 /dr 2 + r 2 + Ar~ 5 / 2 for A > , where the diagonalization 
of H , Eq.(3.3), was carried out in variational spaces of different dimensions n. In Table 3, we 
report our variational computation for upper bounds to the ground-state energy of the Hamiltonian 
H = —d 2 /dr 2 + r 2 + Ar~ 4 along with the eigenvalues reported in the literature. In Table 4, we 
extended our variational analysis to study the Hamiltonian H = —d 2 /dr 2 + r 2 + 1(1 + l)/r 2 + X/r 4 
for A< 1 and for several values of I . We compare our results with the those in the literature, along 
with 'exact' eigenvalues obtained by direct numerical computation of the corresponding Schrodinger 
equation. In order to keep the number of tables of results to a minimum, we first mention the case 
of V(r) = r 2 + which yields the exact energy E = 4 : by using our variational approach we 

obtain an upper bound of E = 4.0000006 for n = 15 with p = 0.73, t = 7.09, and s = 0.01 . The 
reported results in the tables indicate the general usefulness of matrix elements for the investigation 
of the entire spectrum of the spiked harmonic oscillator Hamiltonian for A > and arbitrary a 
in any dimensions and for any angular momentum number I . It is also clear that we don't need 
a very large basis set to produce accurate bounds. It can be seen from the tables that the rate of 
convergence is fast for moderate values of the coupling constant A , while for very small values of 
the coupling constant the rate of convergence is much slower. In general, however, throughout the 
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whole range of values of A , the result from the introduced basis always gives very reliable upper 
bounds. In summary, the basis provides a simple, uniform, and robust variational method. 

4.2 Anharmonic Singular Hamiltonian 

The anharmonic singular Hamiltonians 

H = - + ~ + ar 2 + A + JL , a > 0, c > 0, I = 0, 1, 2, . . . (4.12) 

have attracted considerable attention in part because conditionally exact solutions are possible. 
From the mathematical point of view, this Hamiltonian is a non-trivial generalization of the spiked 
harmonic oscillator (4.3). Znojil [7-8] employed a Laurent series ansatz for the eigenfunctions to 
convert Schrodinger's equation into a difference equation and then used continued fraction solutions 
to obtain exact solutions for the ground-state and the first excited state. Kaushal and Parashar [68] 
simplified Znojil's ansatz to obtain exact ground-state expression 

E = ^fa ^4 + subject to the constraint + b) 2 = c(2l + l) 2 + Sc<Jac. (4.13) 

Guardiola and Ros [37] then used a much simplifier trial wavefunction ip(r) — r^ h l \/c+3)/2 C xp(— r 2 /2— 
v / c/(2r 2 )) for the case of a = 1 and I = to obtain the exact solution for the ground state as 

E — 4:+-^= subject to the constraint condition + b) 2 = c + 8c\/c. (4-14) 

V c 

For example with b = c = 1 , the ground-state is E = 5 and for b = c = 9 , E = 7 , etc. 
Soon afterwards, Landtman [49] performed an accurate numerical calculation and showed that for 
the parameters chosen by Kaushal and Parashar, although the ground-state energy they obtained 
agreed with the numerical calculation, their first-excited energy did not. Varshni [51], in an attempt 
to resolve this problem, obtained four sets of solutions, including one constraint equation for each set 
and showed that the analytic expression for the energy agrees with the numerical result for any one 
among the ground, the first and the second excited states, depending on the particular constraint 
condition satisfied. For higher dimensions, by making use of certain ansatze for the eigenfunction, 
Shi-Hai Dong and Zhong-Qi Ma [52] obtained exact closed-form solutions of (4.12) in two dimensions, 
where the parameters of the potentials a, b, and c again satisfy certain constraints. 

In order to compare our variational results with the exact eigenvalues, we have found for the 
exact ground-state eigenvalue Eq = 5 of the Hamiltonian (4.12) with I = 0,a = b = c = 1 , an 
upper bound of E = 5.000 006 obtained by the diagonalization of a 14 x 14 -matrix. Further, 
the exact energies of 7,7,11,11 corresponding to (a,b,c) — (1, 9, 9), (1, — 7, 49), (1, 45, 225) , and 
(1, —24.5125,600.8623) respectively follow by the optimization of the matrix eigenvalues with ini- 
tial guesses for the variational parameters and matrix dimensions given respectively by (p, t, s) = 
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(1.12, 16.09, 0.11), (1.03, 27.47, 0.07), (1.10, 31. 00, 0.10), and (0.77,40.82,0.01), and 14x14 11x11, 
8x8, and 7x7. These results indicate the generality and the efficiency of our approach. Note 
(b,c) = (-7,49), and (6, c) = (-24.5125,600.8623) also shows the applicability of the method 
in the case of b negative. We further illustrate the applicability of the matrix elements to obtain 
accurate upper bounds to the ground-state of (4.12) for several values of a, b and c. Indeed, for 
(a, b, c) = (1,10,1), (1,10, 10), and (1, 100, 100) , we obtain 6.679 053, 7.138 261 , and 11.791 771 
respectively which results are in excellent agreement with the exact eigenvalues obtained by direct 
numerical integration of Schrodingcr's equation. The precision of the upper bounds to any number 
of decimal places can be achieved by increasing n , the dimension of the matrix. The energies of the 
excited states in arbitrary spatial dimension d are similarly straightforward to find. 



4.3 Perturbed Coulomb Potentials 

Hautot [60], in his solutions of Dirac's equation in the presence of a magnetic field, introduced some 
interesting methods of solving certain second-order differential equations. One of these methods 
deals with the potential operator 

V(r) = -D/r + Br + Ar 2 , A^0. (4.15) 

Hautot obtained exact solutions for only certain relations between the constants A, B, and D . He 
achieved his results by applying the kinetic energy operator to an appropriate wavefunction and 
using the standard procedure of comparing coefficients in the induced recurrence relations. More 
precisely, by introducing [39] 



V>(r) 



exp ( - \ {^Ar 2 + ^0) P^ 1 ' n = Q ^--- 



into the radial Schrodinger equation 

d 2 2d 1(1 + 1) „ D n „ ,\ , . . n 

-TT + --; K —^- L + E+ Br- Ar 2 Wr =0 

dr 1 r dr r z r ) 

One obtains the following three-term recursion relation between the coefficients for (k = 

0,1,2,...): 



(k + 2)(k + 2l + 3) 



D-^={k + 2 + l) 
VA 



E - s/A{2k + 21 + 3) 



AA 



a k = 



This recurrence relation terminates if a k +i = , that is to say E = E„i = \[A{2n + 21 + 3) — |^ 
provided that the parameters A, B, and D satisfies the (n + 1) x (n + 1) -determinant 



dct 



fa b 
ci di bi 

C2 a 2 b 2 



\ 



C n -1 CLn-l bn-1 



[ a k = D- -*L(fc + Z + l), 
0. u-lioiv t b k = (k + l)(k + 2l + 2), 

c fe = E nl - VA{2k + 21 + 1) + ^. 
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For example, we have for the ground-state energy (i.e. k, I = ) of 



V(r) = --+Ar + (Ar) 2 , E = (3 + 2l)A 



4 



(4.16) 



with the ground-state wavefunction given explicitly as ip(r) — r l exp(— i(r + Ar 2 )) . This particular 
case was studied by Killingbeck [17] who obtained the exact solution for the ground-state for j3 > . 
In order to test the variational approach discussed in section 2 and 3, we have employed the matrix 
elements to obtain upper bounds to the exactly solvable cases such as (3 = 0.1, 1, and 2 , we found 
that the upper bounds yields 0.05, 2.75, and 5.75 which are in excellent agrement with the exact 
eigenvalues as given by (4.16). An important consequence of our variational approach are the upper 
bounds that are easily obtained for unconstrained values of D 7 B , and A . In Table 5, we have 
reported our variational results for D = 1 and several values of B and A where we compare 
our results with the upper bounds obtained by the direct numerical integration of Schrodinger's 
equation [39]. In arbitrary dimensions, the matrix elements discussed in Sections 2 and 3 provide 
a uniformly simple, straightforward, and efficient way of obtaining accurate energy bounds for the 
entire spectrum. In order to compare our results with those in the literature, we consider in Table 
6 the radial Schrodinger equation in d -dimensions in the form 



where A = (d + 2£ — 3)/2 . The overall factor of i in the kinetic-energy was incorporate in our 
calculations by multiplying the kinetic energy matrix elements (2.11) by this quantity. To analyze 
the precision of the method, we again compare our results in Table 6 with some special cases for 
which the eigenvalues are known [30] . Results for the excited states within each angular momentum 
subspace (labelled by £) are automatically provided for (up to the dimension of the matrix used), 
and arbitrary spatial dimension dimension d is allowed for in the general expressions for the matrix 
elements. 

4.4. The quartic double-well potential V(r) = —jr 2 + r 4 , 7 > 

The quartic double-well potential 



has a long history of numerical studies (see, for example, [61] and [62] and the references therein). 
Apart from its intrinsic interest, the double- well potential also plays an important role in the quantum 
study of the tunnelling time problem [69], in spectra of molecules such as ammonia and hydrogen- 
bonded solids [70]. Broges et al [71], using supersymmetry techniques, constructed trial wave func- 
tions for variational calculations of the ground-state, first, second, and third excited-states. In their 




(4.17) 



V(r) = -jr 2 + r 4 , 7 > 0, 



(4.18) 
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comparison with the literature, they have used the results obtained from direct numerical integra- 
tion of the corresponding Schrodingcr equation, as reported in [72]. Unfortunately, these numerical 
eigenvalues were not very accurate and the errors are higher than appear in their reported tables. In 
Table 7, we compare our results for the first- and third-exited states with those of Broges et al [71], 
who considered the problem in one dimension; we also include accurate numerical values. 

5. Conclusion 

We have found matrix elements for Schrodinger operators in d spatial dimensions with 
spherically-symmetric potentials of the form V(r) = J2 q a (l) r9 ■ The matrix elements for a given 
angular momentum I are calculated with respect to a finite basis {4>i}™=Q comprising polynomials 
in r with an overall factor of the form r 1 +( t ~ d )/ 2 e - rP . With the inclusion of a scale parameter s , 
the upper estimates are the eigenvalues f]"' (p, t, s) of an n x n matrix eigen equation of the form 
Hv = XNv, where N — [(</>i,</>j)]- For best results, these estimates are optimized with respect to 
the three parameters {p, t, s} for a given n. For the class of problems considered, the basis has 
the advantage that explicit analytic expressions in terms of the Gamma function are available for 
all the matrix elements. The method is robust and flexible enough to yield excellent results for the 
whole class of problems without the need to work with very large matrices. 
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Table 1. Upper bounds E u to the ground-state eigenvalues of H = —-§pz + r 2 + for different 
values of A , a . The eigenvalues for the case a = 2 are obtain for n = 1 and can be compared 
with the exact formula 2 + \/l + 4A . The exponent refers to the dimension (n) of the matrix used 
for the variational computations; the triples in parentheses refer to the approximate initial values of 
the parameters {p, t, s} . The small letters indicate references where the same values were obtained 
in the literature. 



A 


a = 0.5 


a = 1 


a = 1.5 


a = 2 


a = 2.5 


0.0001 


3.000 102 ^ 
(2.00,1.00,1.00) 


3.000 112 ^ 
(2.00,1.00,1.00) 


3.000 138 ^'^ 
(2.00,1.00,1.00) 


3.000 199 980 
(2.00,1.00,1.00) 


3.000 408( 14 < a ) 
(2.00,1.00,0.79) 


0.001 


3.001 022 ^ 
(2.00,1.00,1.00) 


3.001 128 {1 ' b ' d '> 
(2.00,1.00,1.00) 


3.001 382 W 
(2.00, 1.00, 1.00) 


3.001 998 004 
(2.00,1.00,1.00) 


3.004 022 ( 14 < c < e ) 
(2.01,1.06,0.78) 


0.01 


3.010 226 W 
(2.00,1.00,1.00) 


3.011 276 {1 ' b ' d ) 
(1.99,1.00,0.99) 


3.013 794 ( 3 ) 
(1.99,1.00,0.99) 


3.019 803 903 
(1.99,1.01,0.99) 


3.036 744( 15 < c < e ) 
(0.75,1.56,0.01) 


0.1 


3.102 139 ^ 
(2.00,1.00,1.00) 


3.112 067 ^' b -^ 
(1.90,1.00,0.96) 


3.135 053 W 
(1.59,1.00,0.47) 


3.183 215 957 
(2.00,1.18,1.00) 


3.266 874( 18 < c < e ) 
(0.57,1.00,0.00) 


1 


3.009 204 < 13 ) 
(2.18,1.00,0.98) 


4.057 877( 14 < M > 
(0.99,1.00,0.10) 


4.141 893 ( 14 ) 
(1.80,1.06,0.56) 


4.236 067 978 
(2.00,2.23,1.00) 


4.317 311 ( 16 < c < e ) 
(0.69,1.09,0.009) 


10 


12.093 130 (14) 
(1.99,1.00,0.9) 


10.577 483 (14 < b < d) 
(1.00,1.004,0.09) 


9.324 173 (14) 
(2.12,1.69,1.00) 


8.403 124 237 
(1.99,6.40,0.99) 


7.735 111 (6 ' c ' e) 
(1.70,5.85,0.73) 



Ref. [21]. b Ref. [33]. c Ref. [31]. d Ref. [44]. e Refs. [41], [59], and [29]. 



Table 2. Comparison of upper bounds for the ground state energy of the Hamiltonian H = 
— -jpi + r 2 + -^rj2 by different variational techniques. The upper bounds E u are those obtained by 
the present work. The exponent (n) refer to the dimensions of the matrix used for the variational 
computations. 



A 


Ref. [31] 


Ref. [41] 


Ref. [41] 


Ref. [38] 


E u 


0.001 


3.004 075^ 


3.004 074 ( 30 > 


3.004 047 ^ 


3.004 04 


3.004 022 ( 14 > 


0.01 


3.039 409 ( 3 °) 


3.039 244 ( 3 °) 


3.037 474 ^ 


3.037 43 


3.036 744 {15 > 


0.1 


3.302 485 ( 3 °) 


3.296 024 ( 3 °) 


3.269 700 ^ 


3.269 28 


3.266 874 ( 18 ^ 


1 


4.329 449 (30) 


4.32 3 2 63 (30) 


4.318 963 (5) 


4.318 54 


4.317 311 (16) 


10 


7.735 1 36 (30) 


7.735 114^ 


7.735 596 (5) 


7.735 32 


7.735 111 (8) 


100 


17.541 890 ( 3 °) 


17.541 890 ( 3 °) 


17.542 040 ^ 


17.541 92 


17.541 890 ( n > 


1000 


44.955 485 ( 30 ^ 


44.955 485 ( 30 ^ 


44.955 517 ^ 


44.955 49 


44.955 485 W 
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Table 3. Upper bounds for the ground state energy E u for the Hamiltonain H — — A + r 2 + ^ 
for several values of A . We compare our upper bounds E u with the literature. The superscript 
numbers are the dimension of the matrix used; the triples in parentheses refer to the approximate 
initial values of the parameters {p, t, s} . 



A 


E u 


E 


0.0001 


3.022 275 22 
(0.30,3.89,0.00) 


3.022 275 fe 


0.001 


3.068 763 20 
(0.33,8.00,0.00) 


3.068 763°< c , 3.06877° 


0.005 


3.148 352 20 
(0.44,2.36,0.001) 


3.148352 c , 3.14839 h , 3.14835 9 
3.14664°, 3.148352°, 3.05319™ 
3.14835° 


0.01 


3.205 069 20 
(0.44,2.10,0.00) 


3.205067°, 3.20508°, 3.20442 d 
3.20507 fe ,3.20527'\ 3.20507° 
3.205067 d , 3.07522™, 3.23775 fe 
3.20548™, 3.20507°, 3.24980 p 


0.1 


3.575 559 14 
(0.42,9.30,0.00) 


3.575552°< c , 3.57557° 
3.57555 fc , 3.62644 fc , 3.60044^ 


0.4 


4.031 971 22 
(0.50,3.89,0.00) 


4.031971 J ,4.031971 / 


1 


4.494 179 11 
(0.70,11.00,0.00) 


4.494178°, 4.49418 6 
4.49418' c ,4.54879' £ 


10 


6.606 625 14 
(0.44,18.00,0.00) 


6.606623°' c , 6.60662 6 , 6.60662* 
6.64978 fe , 6.609 66 p 


100 


11.265 080 7 
(0.49,72.00,0.00) 


11.265080°, 11.26508°, 11.26508 fe 
11.265 86 p 


1000 


21.369 464 6 
(0.62,100.80,0.00) 


21.369463°'^, 21.36946° 
21.370 26 p 



° Ref. [28]. 6 Ref. [34]. c Ref. [28]. d Ref. [26]. e Ref. [21]. f Ref. [29]. 
s Ref. [20]. h Ref. [2]. 1 Ref. [8]. i Ref. [53]. k Ref. [29]. ' Ref. [20]. 
m Ref. [2]. ™ Ref. [8]. ° Ref. [44]. p Ref. [1]. ° Ref. [25]. s Ref. [38]. 
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Table 4. A comparison between the upper bounds for the Hamiltonian H = —-§pz + r 2 + ^ + ^ , 
for a wide range of values of A = 1(1 + 1) and A , using the present work E u and the bounds 
obtained by Aguilera-Navarro et al [36] (see also [48] and [28]). Accurate numerical results 
E obtained by direct numerical solution of Schodingcr's equation are also presented. The exponent 
refers to the dimension (n) of the matrix used for the variational computations; the triples in 
parentheses refer to the approximate initial values of the parameters {p, t, s} . 



A 


I 


E u a 


E u 


E 


0.001 


3 


9.000 114 279 


9.000 114 279< n ) 


(1.99,7.00,0.99) 


9.000 114 279 




4 


11.000 063 490 


11.000 063 490( n ) 


(1.90,9.00,1.03) 


11.000 063 490 




5 


13.000 040 403 


13.000 040 403< n ) 


(1.98,9.00,0.96) 


13.000 040 403 


0.01 


3 


9.001 142 268 


9.001 142 199< n ) 


(2.10,7.22,0.98) 


9.001 142 199 




4 


11.000 634 795 


11.000 634 788< n ) 


(2.04,6.91,1.00) 


11.000 634 788 




5 


13.000 404 001 


13.000 404 000( n ) 


(2.01,7.00,1.07) 


13.000 404 000 


0.1 


3 


9.011 370 328 


9.011 364 024( 13 ) 


(2.00,5.00,1.19) 


9.011 364 024* 




4 


11.006 336 739 


11.006 336 013< 13 ) 


(2.00,3.99,0.84) 


11.006 336 013* 




5 


13.004 036 546 


13.004 036 433< 8 > 


(1.85,5.97,0.79) 


13. .004 036 433 


1 


3 


9.109 013 250 38 


9.108 657 991 (14) 


(1.70,4.00,0.61) 


9.108 657 991* 




4 


11.062 293 143 4 


11.062 241 722< n ) 


(1.90,8.99,0.78) 


11.062 241 719* 




5 


13.040 025 483 8 


13.040 015 183® 


(1.81,7.19,0.77) 


13.040 015 183 



Table 5. Upper bounds for the Hamiltonian H = — + Br + Ar 2 for different values 

of the parameters B and A . The numerical results in the brackets are the exact eigenvalues as 
obtained by direct numerical integration of Schrodinger equation. The triples in parentheses refer to 
the approximate initial values of the parameters {p, t, s} . 



D 


B 


A 


E u 




1 


2 


3.656 525 (3.657) 
8 x 8,(1.99,1.00,0.73) 




0.1 


1 


1.885 424 (1.885) 
11 x 11,(1.92,1.00,0.75) 




0.5 


1 


2.277 581 (2.278) 
10 x 10,(2.04,1.00,0.86) 




0.1 


0.1 


0.378 305 (0.378) 
12 x 12,(2.21,1.00,1.63) 




0.01 


1 


1.795 268 (1.795) 
8 x 8,(1.08,1.00,0.09) 




0.001 


1 


1.786 212 (1.786) 
8 x 8,(2.04,1.00,0.93) 
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Table 6. Comparison of the eigenvalues for H = —-§pi — — + Br + Ar 2 for different values 
of D , B, and A where E N is calculated from the shifted 1/N expansion [30], the exact su- 
persymmetric values E s [30] and the upper bounds E u obtained by the method of the present 
paper (diagonalization of the dx d matrix elements then minimizing with repect to the parameters 
{p,t,s} ). 



1 


D 


B 


A 


E N 


E s 


E u 





1 


0.447 21 


0.1 


0.171 66 


0.170 82 


0.170 82 6 


1 


1 


0.223 61 


0.1 


0.993 37 


0.993 03 


0.993 04* 


2 


1 


0.149 07 


0.1 


1.509 79 


1.509 69 


1.509 69 3 


3 


1 


0.111 80 


0.1 


1.981 24 


1.981 21 


1.981 21 a 





1 


1.414 21 


1.0 


1.627 56 


1.621 32 


1.621 32 4 


1 


1 


0.707 11 


1.0 


3.411 41 


3.410 53 


4.410 54 s 


2 


1 


0.471 40 


1.0 


4.894 40 


4.894 19 


4.894 19 4 


3 


1 


0.353 55 


1.0 


6.332 78 


6.332 71 


6.332 71 1 





1 


4.47214 


10 


6.226 80 


6.208 20 


6.208 22 4 


1 


1 


2.23607 


10 


11.057 19 


11.055 34 


11.055 34^ 


2 


1 


1.49071 


10 


15.59732 


15.59692 


15.59692 4 


3 


1 


1.11803 


10 


20.093 49 


20.093 36 


20.093 37 iU 





1 


14.14214 


100 


20.753 21 


20.713 20 


20.713 20 4 


1 


1 


7.07107 


100 


35.233 90 


35.230 34 


35.230 34 b 


2 


1 


4.71405 


100 


49.442 67 


49.441 92 


49.441 92 9 


3 


1 


3.535 53 


100 


63.608 60 


63.608 36 


63.608 37 8 









1 


44.721 36 


1000 


66.65904 


66.58204 


66.582 04 3 


1 


1 


22.360 68 


1000 


111.685 01 


111.678 40 


111.678 40 7 


2 


1 


14.907 12 


1000 


156.470 58 


156.469 20 


156.469 20 7 


3 


1 


11.180 34 


1000 


201.215 30 


201.214 87 


201.214 87 7 



Table 7. Upper bounds for the Hamiltonian H = —-^pi — jr 2 + r 4 with different values of 7. 
E\{y) and E$(V) represent the values obtained from the variational method discussed by Broges 
ct al, and E^ and E% are from the present work (with a 10 x 10 -matrix). We have also included 
accurate numerical results E^ and E% obtained by direct numerical integration of Schrodinger's 
equation. 



7 


Ei(V) 


E? 


E» 


E 3 (V) 






0.1 


3.710 64 


3.708 93 


3.708 93 


11.542 58 


11.488 48 


11.488 48 


0.2 


3.618 90 


3.617 01 


3.617 01 


11.386 92 


11.331 27 


11.331 27 


0.3 


3.525 96 


3.523 87 


3.523 87 


11.230 45 


11.173 10 


11.173 10 


0.4 


3.431 79 


3.429 47 


3.429 47 


11.073 07 


11.013 97 


11.013 97 


0.5 


3.336 36 


3.333 78 


3.333 78 


10.914 77 


10.853 87 


10.853 87 


0.6 


3.239 62 


3.236 76 


3.236 76 


10.755 56 


10.692 80 


10.692 80 


0.7 


3.141 55 


3.138 37 


3.138 37 


11.595 47 


10.530 74 


10.530 74 


0.8 


3.042 10 


3.038 56 


3.038 56 


10.434 48 


10.367 70 


10.367 70 


0.9 


2.941 23 


2.937 30 


2.937 30 


10.272 58 


10.203 67 


10.203 67 


1.0 


2.838 91 


2.834 54 


2.834 54 


10.109 78 


10.038 65 


10.038 65 


2.0 


1.726 29 


1.713 03 


1.713 03 


8.433 95 


8.332 87 


8.332 87 



